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ABSTRACT 

If dark matter is composed of neutralinos, one of the most exciting prospects for its 
detection lies in observations of the gamma-ray radiation created in pair annihilations 
between neutralinos, a process that may contribute significantly to the extragalactic 
gamma-ray background (EGB) radiation. We here use the high-resolution Millennium- 
II simulation of cosmic structure formation to produce the first full-sky maps of the 
expected radiation coming from extragalactic dark matter structures. Our map making 
procedure takes into account the total gamma-ray luminosity from all haloes and 
their subhaloes, and includes corrections for unresolved components of the emission 
as well as an extrapolation to the damping scale limit of neutralinos. Our analysis also 
includes a proper normalization of the signal according to a specific supersymmetric 
model based on minimal supergravity. The new simulated maps allow a study of the 
angular power spectrum of the gamma-ray background from dark matter annihilation, 
which has distinctive features associated with the nature of the annihilation process 
and may be detectable in forthcoming observations by the recently launched FERMI 
satellite. Our results are in broad agreement with analytic models for the gamma- 
ray background, but they also include higher-order correlations not readily accessible 
in analytic calculations and, in addition, provide detailed spectral information for 
each pixel. In particular, we find that difference maps at different energies can reveal 
cosmic large-scale structure at low and intermediate redshifts. If the intrinsic emission 
spectrum is characterized by an emission peak, cosmological tomography with gamma 
ray annihilation radiation is in principle possible. 
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INTRODUCTION 

The nature of dark matter in the Universe is one of the most 
fascinating and important mysteries of astrophysics today. 
Based on the evidence gathered during the last decades, we 
know that a suitable particle physics candidate has to be 
neutral and stable (or with a lifetime comparable to the age 
of the Universe) to be the major component of dark matter, 
and that its interactions with ordinary matter have to be 
very weak in order to be consistent with constraints from 
cosmic structure formation. 

Among the different theories that propose a solution to 
the dark matter problem, supersymmetric (SUSY) exten- 
sions of the standard particle physics model with R-parity 
conservation are particularly attractive (see for example 
lMartirj[l998l . for a pedagogical introduction to SUSY). They 
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generically predict the existence of a lightest supersymmet- 
ric particle (LSP), which appears as a plausible dark matter 
candidate. Among the possible LSPs resulting from SUSY 
theories, two are strongly motivated and have been studied 
in detail: the gravitino, which is the spin-3/2 superpartner of 
the graviton, and the lightest neutralino, the spin-1/2 Majo- 
rana fermion that appears in the Minimal Supersymmetric 
Standard Model (MSSM) as the lightest mass eigenstate of 
a mixture of the superpartners of the neutral gauge bosons 
(the fermions B° and W°) and the neutral SUSY Higgs (H° 
and H°). Although gravitinos are interesting viable candi- 
dates t o be the ma j or component of dark matter (see for in- 
stance iFenel 120051 ): Stcffcn (2006) for recent reviews), they 
can be classi fied as extremely weakly interacting particles 
(|Steffenll2007l ) and are therefore very hard to detect exper- 
imentally. On the other hand, the character of neutralinos 
as Majorana fermions (they are their own antiparticle) and 
their somewhat stronger interactions as just weakly interact- 
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ing massive particles (WIMPs) make them more promising 
candidates for detection, possibly in the near future. 

Being WIMPs with expected masses of the order 
of 100 GeV, neutralinos have all the properties required 
to be the cold dark matter (CDM) in the prevailing 
ACDM cosmogony, which is at pres ent the most success- 
ful model of structu re formation (e.g. iKomatsu et al.l 12009 ; 
ISpringel et afll2006h . 

Significant efforts are therefore under- 
taken to find experimental proof for the existence of neu- 
tralinos. These experiments are usually classified as direct 
and indirect searches (for recent reviews on direct and indi- 
rect dark matter searches in the Galaxy s ee iBertonel 120071 ; 
lHoopenl2007l ; ISpooneri [20071 : Ide Boer|[20ol ). 

Direct dark matter search experiments look for the re- 
coil of ordinary matter by scattering of WIMPs in labo- 
ratories on Earth. Some of the experiments that have been 
looking for signals are: CDMS (lAkerib et al.ll2004h , XENON 
I Angle et al.1l200Sft. ZEPLIN (|Alner et alj|2005l ) and DAMA 
i Bernabei et al.ll2000l . I2008T ). The latter is the only experi- 
ment that has reported a positive signal so far, a finding that 
has not yet been confirmed by another experiment, however. 

The interaction rate of neutralinos with ordinary nu- 
clei fundamentally depends on two factors, the flux of neu- 
tralinos on Earth and the interaction cross section of neu- 
tralinos with a given nuclei. The latter can be computed 
theoretically within the context of the SUSY model, but 
the former depends strongly on the phase space distribu- 
tion of dark matter at the Solar circle, which needs to be 
predicted by cosmological models of structure formation. 
In standard analysis, a constant dark matter density and 
a simple Maxwell-Boltzmann distribution for the particle 
velocities is often assumed, properly modified to include the 
motion of the S un around the Gala xy and that of the Earth 
around the Sun (jFreese et al.ll200"lh . However, it is not at all 
clear that this simple assumption is justified, because cold 
dark matter is expected to produce a large number of over- 
lapping cold streams locally, and may leave residual struc- 
ture in energy space. High-resolution N-body simulations 
of structure formation combined with an explicit modelling 
of the local stream density evolution provide one promising 
new approach to pro vide a more accurate characte rization 
of the possible signal (|Vogelsberger et"al . 20081. |2009\ Other 
dynamical local Solar Sys tem effects may also be important 
l|Lundberg fc Edsio|[20oi ). 

Indirect dark matter searches look for the products 
of WIMP annihilations generated in cosmic structures. 
Such poss ible annihilation pr oducts incl ude, for example, 
posit rons (|Baltz fc E dsio 1999), neutrinos |Berezinskv et all 
1996) and gamma-ray photons. Interestingly, the satel- 
lite PAMELA has recently reported an anomalously high 
positron abundance in cosmic rays at high energies 
l)Adriani et al.l I2009T ) . Although this signal could be ex- 
plained by a number of astrophysical sources, dark matter 
annihilations offer a particularly attractive explanation un- 
der certain scenarios (jBergstrom et ai]l2008h . and this possi- 
bility has created a flurry of activity in the field. In a related 
experiment, the FERMI satellite has measured an excess 
in the total flux o f electrons and positrons in cosmic rays 
l|Abdo et alj|2009\ which, although is smaller than those 
reported by the ba l loon experiments ATIC and PPB-BETS 
ijChang et al.ll200a : iTorii et al.ll2008l ) is nevertheless signifi- 
cant. 



Within our Galaxy, the annihilation signal is ex- 
pected to be particularly strong in the high density re- 
gions around the galacti c centre (|Bertone fc MerrittJ 12003 ; 
Ijacholkowska et al. 2006). However, since other astrophys- 
ical sources contribute prominently to the gamma-ray flux 
there, it may be difficult to disentangle a possible neutralino 
annihilation signal from the rest of "ordinary" sources. 
A more promising approach may lie in an analysis of 
the gamma-ray radiati on produced in th e smooth Milky 
Way halo as a wh ole dStoehr et alJ 120031: I Fornengo et al.l 
12004 Ide Boer et all |2004 : ISpringel et al.l |200ST ). Particular 
attention has al so been given to sub s tructures wi thin the 
MW halo (e.g. ISiegal-Gaskinsl 120081 : lAndol 120091 ). which 
have been suggeste d to dominate the ann ihilation signal 
observed at Earth ijBerezinskv et all l2003h , even though 
they have recently been shown to be more difficult to de- 
tect compared with the signal from the galactic centre 
jSpringel et alj2008h . Searches in other nearb y galaxies have 
also been considered, su ch as Andromeda (|Falvard et al.l 
|2004 lLavalle et all 120061) and other galaxies in the Local 
Group ( Wood et alj|2008h . Clus ters of galaxies may also be 
interesting targets for searches l|Jeltema et al.ll2009f ) due to 
their large content of dark matter, even though a variety of 
other possible sources of high energy gamma-rays in clus- 
ters, such as cosmic rays generated in structure formation 
shocks, active galactic nuclei and supernovae, can lead to 
confusion and may make it challenging to cleanly separate 
the dark matter signal from the confusing sources. 

On still larger scales, we expect that traces of dark 
matter annihilation should be present in the extragalac- 
tic gamma-ray background (EGB) radiation, which con- 
tains gamma-rays coming from all possible cosmic sources. 
In the present work, we will focus on these extragalactic 
gamma-ray photons. The EGB has been measured by differ- 
ent gamma-ray satellites, in particular in the energy range 
between 1 MeV an d 30 GeV by COMPTEL and EGRET 
|Strong et al. 1 120041 ). In the near future, our understanding 
of the EGB will greatly improve by the Large Area Tele- 
scope aboard the recent ly launched F ERMI satellite (for- 
merly known as GLAST. iAtwood et al.ll2009h . which covers 
an energy range between 20 MeV and 300 GeV and features 
a greatly improved sensitivity compared with its predeces- 
sor EGRET. Several sources of different astrophysical origin 
are known to contribute significantly to the EGB, includ - 
ing blazars dAndo et al.ll2007al;lKneiske fc Mannheimfe oOS). 
type la supernovae (|Strigari et al.ll2005l 1. cosmic ra ys accel- 
erated at cosmological structure formation sh ocks (|Miniatil 
120021 : iMiniati et alj|2007l:|jubelgas et alj|2008h . among oth- 
ers (see iDermerl 120071 ; iKneiskel 120081 . for recent reviews). 
Of course, dark matter annihilation also contributes to the 
EGB, and it is in this component that we focus on in this 
work. 

Using standard analytic approaches to describe cos- 
mic structure formation, like halo mass functions and clus- 
tering amplitude as a function of redshift, it is possible 
to make predictions for the relative gamma-ray flux com- 
ing from halo es and their subhaloes as a function of mass 
and redshift (|Ullio et al.ll2002l : iTavlor fc Silklbool ) and to 
use this information to estimate the contribution of dark 
matter annihilation to the EGB. Both, the isotropic and 
anisotropic components of this rad iation have been pre- 
dicted in this way in previous works |Elsasser fc Mannheim! 
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20051; lAndo fc Komatsul 120061; ICuoco et all 120071 , 120081 ; 
de Boer et all 120071; iFornasa et al l 120091 ). In particular, 
Ando fc Komatsu (2006) first pointed out that the density 
squared dependence of dark matter annihilation should leave 
a specific signature in t he angular power spectrum of the 
EGB. We also note that ISiegal-Gaskins fc Pavlidoul (|2009l ') 
have shown that the energy dependence of the angular power 
spectrum of the total diffuse gamma-ray emission can be 
used to disentangle different source populations, including 
that associated with dark matter annihilation. 

However, it is not clear how accurately these analytic 
approaches capture the non-linear structures resolved in the 
newest generation of high-resolution cosmological simula- 
tions. Also, the previous analysis have so far been restricted 
to statistical statements about the power spectrum of the 
EGB, or its mean flux. For a full characterization of the 
signal, it would be useful to have accurate realizations of 
maps of the expected gamma-ray emission over the whole 
sky. Such maps contain in principle the full information that 
can be harvested, and allow, for example, correlation studies 
with foreground large-scale structure in terms of galaxies or 
dark matter. 

In the present paper, we therefore focus on techniques 
to predict the extragalactic gamma-ray background directly 
from cosmological N-body simulations. To this end, we use 
the Millennium-II simulation. With 10 10 particles in a ho- 
mogeneously sampled volume of (100 /i _1 Mpc) 3 , it is one of 
the best resolved structure formation simulations to date. 
We use a ray-tracing technique to accumulate the signal 
from all haloes and subhaloes over the past light-cone of 
a fiducial observer, positioned at a plausible location of the 
Milky Way in the simulation box. Through careful account- 
ing for unresolved structure, we can extrapolate our signal 
to the damping scale limit and obtain the first realization 
of realistic full sky maps of the expected gamma-ray radia- 
tion from dark matter annihilation. The maps can be used 
both to check earlier analytic estimates, and to develop new 
strategies for identification and exploitation of the annihi- 
lation component in the EGB. For example, we show that 
difference maps at different energies could be used to un- 
cover nearby cosmic large-scale structure, which could be 
correlated with other probes of large-scale structure to more 
unambiguously identify the origin of this background com- 
ponent. 

This paper is organized as follows. In Section 2, we 
briefly review the analytic formulae associated with the con- 
tribution to the EGB from dark matter annihilation. In Sec- 
tion 3, we describe the SUSY model chosen for our analysis. 
We then present the N-body simulation used, and the ef- 
fect of dark matter clustering in the annihilation's gamma 
ray flux in Section 4. In Section 5 we describe our method 
to generate simulated sky maps of the EGB radiation from 
dark matter annihilation, and we analyse the energy and an- 
gular power spectra of the signal. We give a summary and 
our conclusions in Section 6. 



2 DARK MATTER ANNIHILATION RATES 
AND THE EGB 

If dark matter is made of particles that can annihilate with 
each other, then the number of annihilations per unit time 



in a given volume V is given by (e.g. IGondolo fc Gelminil 

Eel]): 



2m\ 



I P 2 x {x)d 3 x, 
Jv 



(1) 



where m x and p x are the mass and density of the dark mat- 
ter species, and {av) is the thermally averaged product of 
the annihilation cross section and the M0ller velocity. 

Since the particle physics details related to the intrin- 
sic nature of dark matter are encoded in the term outside 
the integral in Eq. {T}, it is convenient to remove this de- 
pendence in analyzing the effect of dark matter clustering 
on the annihilation rate. For this purpose, it is customary 
to introduce the ratio of gamma ray emission coming from 
a given halo of volume V to the average emission of a ho- 
mogeneous distribution of the same amount of dark matter 
within this volume: 



Ph v Jv 



(2) 



where ph is the average density of dark matter particles in- 
side V . 

For instance if the region of interest is a dark matter 
halo, then Eq. ([2} can be solved analytically by assuming a 
specific density distributi on. In particular, fo r the spherically 
sy mmetric NFW profile (|Navarro et al.f l997'). fh(V) is given 
by l|Tavlor fc Silkll2003l ): 



fh(CA) 



cj [1 - 1/(1 + cj)] 
9[ln(l + c A ) -ca/(1 + ca)] 2 



(3) 



where ca = r&/r s is the concentration parameter, r s is the 
scale radius of the NFW profile, and ta is the limiting radius 
of the halo, connected to its enclosed mass Ma through: 



4 _ 3 4 
Ma = g^p/i^A = g7rAp c 



3 

itr A , 



(4) 



where A is the overdensity parameter. Different values for A 
are used in the literature, with the two main choices being 
A = 200 ijNavarro et al.l 1 1997ft . and A ~ 178f2^ 45 for a 
ACDM model with Q. m + fi A = 1 l|Eke et alj|200ll ). 

For spherically symmetric haloes, t he integral in Eq. Q 
can also be written as the scaling law ijSpringel et al.ll2008T ) 
L' h oc VmaxAhaif, where V max _ is the maximum rotational 
velocity of the halo and rhaif is the radius containing half 
of the gamma-ray flux. For a NFW halo, the latter formula 
reduces to 



r' f 2 / \ 1-23 V x 

L h = / pNFw{r) dV = —2 — 



4 

max 



(5) 



where r max is the radius where the rotation curve reaches 
its maximum. We will use Eq. ([5]) for our analysis in this 
paper. The most recent high-resolution N-body simulations 
support a slightly revised density p rofile that becomes grad- 
ually shallower towards the centre ((Navarro et al.|[2009h , de- 
viating from the asymptotic power-law cusp of the NFW 
profile. For this Einasto profile, the coefficient 1.23 in Eq. ([5]) 
changes to 1.87, an insignificant difference compared to other 
uncertainties in the analysis of the dark matter annihilation 
radiation. 

We can extend Eq. ((2]) to larger volumes containing sev- 
eral haloes of different sizes. If we assume that all dark mat- 
ter particles in a given cosmic volume Vb belong to virialized 
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haloes, then the ratio of gamma-ray emission coming from 
all haloes contained in Vb with masses larger than M m i n 
to the emission produced by a smooth homogeneous distri- 
bution of dark matter, with average density pB, filling this 
volume is given by: 



f(M > M min ) = 



1 



PbVb 



ANv B (M)plV{M)f h (y{M)) 
(6) 

where ANv B (M) is the number of haloes in the volume Vb 
with masses in the range [M, M + dM], which can be written 



dNvAM) = Vb^AM = fe^gdlogM, (7) 



dlogM 



with dn(A/)/dlogA/ denoting the differential halo mass func- 
tion. Thus: 



f(M > M, 
- f 



(dloiS) plV(M)f h (V(M))d\ogM. 



(8) 



In the literature, f(M > 
Hie "flux multiplier" (Eq. 18 o f iTavlor fc Silk! 120031) o r 



Af m j n ) has been ca lled 



the "clumping factor" (Eq. 6 of lAndo fc Komatsul 12006). 
The limiting lower mass in Eq. (|8]) is the minimum mass 
scale for virialized haloes surviving today, which is roughly 
given by the free streaming length of dark matter parti- 
cles: ~ lQ- e Mq f or m x = 100 GeV ijHofmann et all (200 it 
iGreen et al.ll20oj ). see section 4.1 for a discussion on this 
cutoff mass. Although many of the small scale haloes have 
suffered from tidal truncation or were disrupted, it is pos- 
sible that they contribute significantly to the gamma-ray 
emission produced by dark matter annihilation. 

Eqs. ([2]) and ([8} allow us to study the gamma ray flux 
coming from individual haloes or groups of haloes in a given 
region of the Universe at a given redshift (both can be ap- 
plied to any redshift, provided that p x , pu and ps are com- 
puted at such redshift). 

Following Eq. we define the gamma ray emissivity 
(energy emitted per unit energy range, unit volume and unit 
time) associated with dark matter annihilation as 



E. 



d7V 7 (av) 
' dE~, 2 



Px 



(9) 



where dN-y/dE^ is the total differential photon spectrum 
summed over all channels of annihilation. Of all the quan- 
tities in this equation, only p x depends on the spatial dis- 
tribution of dark matter, the rest is related to its intrinsic 
properties. For the analysis of the latter, since we consider 
dark matter to be made of neutralinos within a SUSY see- 
conveniently define the so called SUSY factoiQ: 



nario, we < 



fa 



dNj (av) 
dE-y ml 



(10) 



The signal that we are ultimately interested in is the 
contribution of neutralino annihilation to the EGB, or more 
specifically, the spatial distribution of the gamma-ray radi- 
ation coming from dark matter annihilation integrated over 



1 This term has been extensively used previously in the literature, 
for example in iFornengo et al.l l|2004l ) , /sus Y is related to the 
integral quantity <J>SUSY defined there. 



all redshifts along the line-of-sight for a fiducial observer, lo- 
cated at z — 0, for all directions on its two-dimensional full 
sky. The values we want to calculate for each pixel of such 
a sky map are that of the specific intensity, the energy of 
photons received per unit area, time, solid angle and energy 
range: 



I-,,o = J e 7 , (£ 7 ,o(l + z), 



z) 



dr 



(1 + *) 



(11) 



where the integral is over the whole line of sight, r is the co- 
moving distance and E lt o is the energy measured by the ob- 
server at z — 0. Note that e 7 ,o is evaluated at the blueshifted 
energy (1 + z)E 1: o along the line-of-sight to compensate for 
the cosmological redshifting. 



3 DARK MATTER ANNIHILATION AND THE 
SUSY FACTOR 

In the present section we describe the particular SUSY 
model that we used to compute the gamma-ray emissiv- 
ity produced by dark matter annihilation. In particular, we 
provide a closer examination of the SUSY factor given in 
Eq. (HSJ. 

A detailed particle physics model for the nature of the 
neutralino as dark matter, and its mass and interactions, is 
needed to compute /susy- We restrict our analysis to the 
framework of th e minimal supergrayity m odel (mSUGRA, 
for reviews see iNillesI 1 1984 iNathl I2003T ). The mSUGRA 
model is popular, among other reasons, due to its relative 
simplicity, which reduces the large number of parameters in 
the general MSSM (> 100) to effectively four free parame- 
ters, mi/2, mo, tan/3, Ao and the choice of sign for p. The 
parameters m^, mo and Ao are the values of the gaugino 
and scalar masses and the trilinear coupling, all specified at 
the GUT scale, tan/3 is the ratio of the expectation values 
in the vacuum of the two neutral SUSY Higgs (H° y H%), 
and finally, p is the Higgsino mass parameter. 

The general 5-dimensional parameter space of the 
mSUGRA model is significantly constrained by various re- 
quirements: consistency with radiative electroweak symme- 
try breaking (EWSB) and experimental constraints on the 
low energy region, which are given by current lower mass 
bounds for the lightest SUSY Higgs after EWSB, ra/,, as well 
as concordance between experimental and theoretical pre- 
dictions for the decay b — > sy and the anomalous magnetic 
moment of the muon (gr M — 2)/2. Importantly, it is normally 
also assumed that the major component of dark matter is 
the LSP of the theory, then the LSP must be neutral and 
f^LSP = Hdm , which reduces the allowed regions in the pa- 
rameter space considerably. For typical mSUGRA scenarios, 
the LSP that satisfies these conditions is the lightest neu- 
tralino. Its relic density fi x needs to be properly calculated 
by solving the Boltzmann equation, with the parameters of 
the theory chosen such that the desired dark matter density 
results. 

Observational constraints on the abundance of dark 
matter, J7dm, come from various sources, including the cos- 
mic microwave background (CMB), galaxy clusters, super- 
novae, Lyman-a forest data, and weak gravitational lensing. 
For the purposes of the analysis in this section, we take the 
la bound for the amount of dark matter in the Universe 
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Figure 1. Sketch showing the allowed parameter space of the 
mSUGRA model when the constraint on the relic density accord- 
ing to Eq. M2\ (dark blue region) is imposed. The experimental 
constraint on the mass of the lightest SUSY Higgs, mj,, is also 
shown for reference. The grey region in the lower right corner is 
where the LSP is charged, whereas in the grey region in the upper 
left corner there is no radiative EWSB. The abbreviated names 
for various important regions in the parameter space are marked 
in the figure with capital letters, see text for details. 

according to the five year results of WMAP data combined 
with distance measurements from type la superno vae and 
Baryon Acoustic Oscillations (|Komatsu et alj|2009l ): 

0.1109 sC !!dm/i 2 < 0.1177. (12) 

The general picture for the allowed regions in the pa- 
rameter space of the mSUGRA model resulting from the 
different constraints can be qualitatively described in the 2- 
dimensional plane (mo — 7711/2), as shown in Fig. [1] The blue 
area in the figure represents the allowed region according to 
Eq. (|12[) for the observational bounds on the abundance of 
dark matter. 

The different regions in the mSUGRA parameter space 
that satisfy thes e constraints have been studied abundantly 
in the past (e.g. iBaer fc Balazsll2003l ; lBattaglia et al.ll20ol 
iFend 12005; Bcl anger et al.ll2005h and have received generic 
names, marked with capital letters in the figure: (B) bulk 
region (low values of mo and 77I1/2), (FP) focus point re- 
gion (large values of mo), (CA) co- annihilation region (low 
mo and m x < m-f x , where f\ is the lightest slepton) and 
(RAF) rapid annihilation funnel region (for large values of 
tan/3 and a specific condition for mjj). The regions in grey 
(light) represent zones where the TSP is charged (lower right 
corner) and where there is no radiative EWSB solution (up- 
per left corner). The dashed line is shown as a reference for 
the experimental constraints for mj,. 



2 The relation which should hold is: 2m Y £ s m^, where A is the 
CP-odd Higgs boson; see for example iFengl (2005) for details. 



Our goal in this section is to analyze the supersymmet- 
ric factor /susy predicted by the mSUGRA model in the dif- 
ferent allowed regions depicted in Fig. [1] To do so, we define 
several "benchmarks points" which fulfill all the experimen- 
tal constraints discussed above and which can be taken as 
representative of the different allowed regions in parameter 
space that we just described. T he benchmark points were 
selected following the analysis of lBattaglia et all (12004 ) and 
iGondolo et al.l |2004 ) , but slightly modified to be consistent 
with the cosmological c onstraint of Eq. (|12|1. We u se the nu- 
merical code Dark SUSY dGondolo et al.ll2004l20u5l ) with the 
interface ISAJET |Baer et al.ll2003l ) to analyse these bench- 
mark points. Table 1 lists the 11 selected benchmark points 
with their corresponding values for the five mSUGRA pa- 
rameters, the value of the lightest neutralino's mass, its relic 
density, and the name of the region they belong to in Fig. [T] 

Fig. [2] shows the energy spectrum of /susy computed 
for four of the benchmark points presented in Table 1: A 
(black dashed line), E (red dotted line), K (green dash- 
dotted line), and L (blue solid line), which are represen- 
tatives of the regions CA, FP, RAF and B, respectively. 
The main features of Fig. [2] are determined by the photon 
spectrum d-/V 7 /di? 7 which receives contributions from thre e 
mechanisms of photon production l|Bringmann et all feoOS): 
(i) gamma-ray continuum emission following the decay of 
neutral pions produced during the hadronization of the pri- 
mary annihilation products at the tree-level; this mechanism 
is dominant at low and intermediate energies (in Fig. [5] it 
dominates the spectrum in most of the energy range); (ii) 
monoenergetic gamma-ray lines for neutralino annihilation 
in two-body final states containing photons, which is allowed 
in higher order perturbation theory; the most relevant final 
states containing these lines are XX ~~ ^ 77 an d XX ~> ^°7, 
where Z° is the neutral Z boson; for non-relativistic neu- 
tralinos the energy of each outgoing photons in these pro- 
cesses is E-y = m x and E 1 = m x (l — m^o/^m 2 ,), respec- 
tively (in Fig. [2] these lines are prominent for benchmark 
point E (red line) and negligible in the rest of the cases); 
(iii) when the final products of annihilation are charged, in- 
ternal bremsstrahlung (IB) needs to be considered, it leads 
to the emission of an additional photon in the final state; 
the IB contribution to the spectrum is typically dominant 
at high energies resulting in a characteristic bump (bench- 
mark points A and L have an important contribution from 
IB at energies close to m x , for the other two points the 
contribution from IB is negligible). Processes (ii) and (iii) 
are subdominant relative to process (i), but they display 
distinctive spectral features intrinsic to the phenomenon of 
annihilation. 

Additionally, the annihilation of neutralinos into elec- 
trons and positrons contributes indirectly to the high 
energy photon spectrum when low energy background 
photons, such as those in the cosmic microwave back- 
ground (CMB), are scattered to higher energies via inverse 
Compton (IC) scattering by these electrons and p ositrons 
|Profumo fc Jeltemallioog : iBelikov fc Hooperll2009l ). In the 
case of CMB photons, the energy peak of this IC photon 
emission {E^ c ) is approximately independent of redshift 
and is given by: Elf w E^ MB {E e ± /m e ±c 2 ) 2 , where E% MB 
and E e ± are the average energies of the CMB photons 
and e respectively. For neutralino masses in the 100 GeV 
range, this energy peak lies in the X-ray regime (~ 10 keV) 
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Model" ACDEFGH IJKL 



m l/2 


550 


435 


520 


399 


811 


396 


930 


395 


750 


1300 


450 


mo 


114 


94 


113 


2977 


4307 


118.5 


242.5 


193 


299.5 


1195 


299 


tan/3 


7 


11 


10 


30 


30.4 


20 


20 


35 


35 


46 


47 


sign(^) 


+ 


+ 




+ 


+ 


+ 


+ 


+ 


+ 




+ 


m x 


226.4 


175.1 


214.9 


154.2 


337 


160.4 


394.2 


160.6 


315.5 


566.1 


184.9 


Q x h 2 


0.115 


0.114 


0.116 


0.118 


0.112 


0.113 


0.118 


0.118 


0.114 


0.116 


0.111 


Region^ 


CA 


B 


CA 


FP 


FP 


B 


CA 


B 


CA 


RAF 


B 



Table 1. Selected benchmark points in the mSUGRA parameter space. For our main analysis, we will focus on point 'L', which 
predicts an annihilation flux close to the maximum possible for the constrained mSUGRA models. 

a The letters to represent each point were chosen to follow the notation of Battaglia et al. | Hool) and lGondolo et al. | [120041) 
^ B for Bulk region, FP for Focus Point region, CA for Co-annihilation region and RAF for Rapid Annihilation Funnel region 
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10" 36 r % \ \ 

1 cr 37 L i lij . J 

0.1 1.0 10.0 100.0 

E, [GeV] 

Figure 2. Supersymmetric factor as a function of 7-ray emission 
energy for a selection of the benchmark points described in Table 
1: A, E, K and L, which are shown with black (dashed), red 
(dotted), green (dash-dotted) and blue (solid) lines, respectively. 



llProfumo fc Jeltemai r2009). Unless the neutralino mass is 
larger than 1 TeV and with a hard e ± spectrum, the con- 
tribution of this mechanism is not significant in the energy 
range relevant to the present study (> 0.1 GeV). We note 
that we are currently working on the analysis of the X-ray 
extragalactic background radiation from dark matter anni- 
hilation and we will present our results in a future work. 

In Fig. [2] we have shown only 4 of the 11 benchmark 
points we considered; the rest show very similar features 
and lie in between these 4 cases. Fig. [2] also indicates that 
the normalization of the spectrum of /susy typically varies 
by three orders of magnitude among the different regions 
of the allowed parameter space. Also there appears to be a 
trend between the different allowed regions in the mSUGRA 
model: /susy has the largest value for benchmark points 
on the focus point region (FP, red dotted line), followed 
by the bulk region (B, blue solid line), the rapid annihila- 



tion funnel region (RAF, green dash-dotted line), and finally 
the co-annihilation region (CA, black dashed line). This last 
feature could be of importance for constraining the allowed 
mSUGRA parameter space even more: in principle, a pre- 
cise measurement of the gamma-ray flux coming from dark 
matter annihilation could discriminate between the charac- 
teristic regions FP, B, CA and RAF, considering for example 
that the difference in /susy between the CA and FP regions 
is around two orders of magnitude. However, such inferences 
will only become possible if the remaining sources of uncer- 
tainty can be reduced to less than two orders of magnitude. 
This is a demanding goal, given the limited knowledge we 
have on some of the other astrophysical factors that play an 
important role in the production of these gamma-rays, as 
well as the observational difficulties in properly subtracting 
from an observed gamma-ray signal the contributions from 
astrophysical sources unrelated to dark matter annihilation. 
Nevertheless, 'dark matter astronomy' remains an interest- 
ing possibility in the light of the results shown in Fig. [2] 

In our subsequent analysis, we will now choose one par- 
ticular benchmark point and adopt its /susy value for the 
rest of our work, keeping in mind the results above. Our cho- 
sen benchmark point is the model L (that we highlight with 
a thick blue solid line in Fig. [2]), which gives an upper limit 
on /susy for the benchmark points in the bulk region and is 
close to the maximum value of /susy for all the benchmark 
points analyzed. The latter is the main reason to choose 
this particular model because it is desirable that fsuSY, 
and hence the predicted gamma-ray flux, has a large value 
among the different theoretical possibilities. This optimistic 
choice as far as prospects for future detections are concerned 
then also demarcates the boundary where non-detections 
can start to constrain the mSUGRA parameter space. But 
there are other reasons as well. The mass of the neutralino 
in this model is "safely" larger than the lower mass bounds 
coming fro m experimental con straints (m x > 50 GeV ac- 
cording to iHeister et alj 120041 ') , but low enough to be de- 
tectable in the experiments available in the near future. Also, 
a high value of the parameter tan/3 seems to be favor ed by 
other theoretical expectations (e.g. iNrinez et alj|2008h . We 
note that although benchmark point L does not have promi- 
nent gamma-ray lines in its spectrum, it does have an im- 
portant IB contribution and serves the point of showing the 
implications of a peak in the annihilation energy spectrum 
for the EGB. 
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We note that recent analyses on the cosmic-ray 
anomalies as measured by PAMELA and FERMI favour 
higher neutralino masses (~ 1 TeV) when these anoma- 
lies are explained mai nly by dark matter annihilation (e.g. 
iBergstrom et al.ll2009l ). However, it is important to keep in 
mind that other astrophysical sources, such as pulsars and 
supernovae remnants, are expected to have a meaningful 
contribution to the solution of the cosmic ray anomalies. At 
present, these experiments are only suggestive of dark mat- 
ter particles in the TeV range, not definite. Our goal is to 
give an account of the general features of the gamma-ray 
spectrum from annihilation and their imprint on the EGB. 
For this purpose, our fiducial model with m x ~ 200 GeV is 
adequate. Models with higher masses would shift the energy 
scale to higher values, but our general conclusions would 
remain unaltered. 



4 GAMMA-RAY LUMINOSITY FROM 
ANNIHILATION IN HALOES AND 
SUBHALOES: THE ASTROPHYSICAL 
FACTOR 

In this section we analyse the astrophysical part of the anni- 
hilation gamma-ray flux, arising from the clustering of dark 
matter into haloes and subhaloes across cosmic time. To this 
end we analyse the dark matter haloes in one of the most 
recent state-of-the-art cosmological N-body simulations. 

Specific ally, we use the "M il lenniu m-H" simulation 
(MS-II) of iBovlan-Kolchin et alj (120091 ) which has the 
same particle number (2160 3 ) and cosmological parameters 
(Sl m = 0.25, Sl A = 0.75, h = 0.73, cr 8 = 0.9 and n s = 1, 
where Sim and SIa are the contribution from matter and 
cosmological constant to the mass/energy density of the 
Universe, respectively, h is the dimensionless Hubble con- 
stant parameter at redshift zero, n s is the spectral index 
of the primordial power spectrum, and as is the rms am- 
plitude of linear mass fluctuations in 8 h, _1 Mpc s pheres at 
redsh ift zero) as the Millennium simulation (MS-I) (|Springell 
l2005h but with a box size that is 5 times smaller, equal to 
L = 100 /i _1 Mpc on a side, thus having a mass resolution 
of 6.89 x 10 6 /i _1 M Q , 125 times smaller than in MS-I. Typ- 
ical Milky-Way sized haloes are resolved with several 10 5 
particles, while clusters of galaxies have about 50 million 
particles. 

For the dynamical evolution of the MS-II, a fixed co- 
moving gravitational softening length of e = 1 /i _1 kpc 
(Plummer-equivalent) was used. The MS-II was processed 
on-the-fly to find haloes using a friend-of-friends (FOF) al- 
gorithm, followed by a post-processing step to identify gravi- 
tationally bound dark matter substructures down to a limit 
of 20 parti cles within each gro up based on the SUBFIND 
algorithm dSpringel et al.l l200lh . A given FOF halo is de- 
composed by SUBFIND into disjoint subgroups of self-bound 
particles. The most massive of these subgroups represents 
the "ma in" or "background" halo of a given FOF halo (see 
Fig. 3 of lSpringel et al.ll200ll ). For the reminder of this work 
we will refer to these subgroups as main haloes or simply 
haloes, and to the rest of the identified subgroups of a given 
FOF halo as subhaloes of that main halo. SUBFIND also 
computes several properties for each subgroup; the ones we 
will use extensively in this work are Knax and r max , the max- 



imum circular velocity of a (sub)halo and the radius where 
this maximum is attained. We also use ri/%, the radius which 
encloses half of the mass of a subhalo, and T2oo, the radius 
where the mean density of a main halo is 200 times the mean 
background matter density. 

Even with the large particle number available in the 
MS-II simulation, it is important to note that a direct eval- 
uation of halo luminosities by estimating local dark matter 
densities at the positions of each simulation particle (for 
example based on the SPH kernel interpolation technique) 
would be seriously affected by resolution limitations. This 
is because most of the emission originates over a tiny ra- 
dial region, close to the very centre of a halo, where the 
density is easily underestimated even for otherwise well re- 
solved structures. In fact, it requires of order 10 s — 10 9 par- 
ticles to obtain a converged estimate of the emission from 
the smooth par t of a cuspy dark m atter halo in a brute- 
force approach dSpringel et al.ll2008l ). On the other hand, it 
is numerically very much easier to resolve the existence of a 
halo or subhalo, and to accurately estimate its total mass, 
as well as its primary structural properties such as Vmxc 
and r max . Our strategy will therefore be to adopt an ana- 
lytic fit for the density profile of each detected dark matter 
halo and to predict its expected emission by integrating this 
profile. In other words, we will use the scaling relation ((5)1 
to compute L' h for main haloes and subhaloes alike. Since 
Vinax and r max can be measured reasonably accurately even 
for poorly resolved haloes, this gives us a reliable account- 
ing for the total emission of all haloes down to a fairly low 
mass limit. Importantly, we can completely avoid in this way 
to introduce strong mass- dependent resolution effects in our 
measurements of luminosity as a function of halo mass. 

For simplicity, we shall assume that all haloes have a 
NFW structure, where the local logarithmic slope of the 
density profile tends asymptotically to —1 for small radii, 
even though the Aquarius project, a recent set of simula- 
tions aimed at following the formation of MW-sized haloes 
at different mass resolutions, demon strated that the Eina sto 
profile provides an even better fit l|Navarro et alj [20091 ). a 
finding that is also supported b y the independent 'GHALO' 
simulation (|Stadel et al.l 2009). In the Einasto profile, the 
local logarithmic slope of the density profile changes with ra- 
dius according to a power law with exponent a s h (typically 
called shape parameter). According to the Aquarius project, 
a 3 h is of the order .16 — 0.17 an d varies from hal o to halo 
l|Navarro et al.ll2009l ). Previously. [Gap et all (|2008h found a 
dependence on halo mass for the value of the shape param- 
eter (going from 0.16 for galaxy-haloes to 0.3 for massive 
clusters). However, for the purposes of our work, the NFW 
profile is a good enough approximation and using it sim- 
plifies comparisons to a number of results in the literature. 
In any case, we note that using the Einasto profile (with 
a typical shape parameter of 0.17) instead would increase 
our results for the fluxes by ~ 50%. This value changes only 
very slightly over the range of values for a s u found for dark 
matter haloes. Such an uncertainty on the dark matter den- 
sity profile is considerably lower than others in our analysis, 
namely, the photon spectrum from annihilation (see Fig. [2]| 
and the contribution from unresolved haloes and subhaloes 
(see sections 5.3 and 5.4). 

Nevertheless, it is important to take into account the 
impact of finite numerical resolution in the measured val- 
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ues of Vmax and r max for the smallest haloes that c an be re- 
solved. As has been shown bv lSpringel et all (|2008l ) . the val- 
ues of r max are increasingly overestimated for smaller haloes 
whereas for Vrn ax the opposite holds. A partial correction for 
these numerical effects can be obtained by trying to account 
for the effect of the gravitational softening e for haloes with 
values of r max close to e. Assuming adiabatic invariance of 
orbits between softened and unsoftened versions of the same 
halo, one obtains to leading order: 



^"m aX — ?"m aX [l "1" (t/fmax) ] 
V" ma x = VUw[l + (<0-ma X ) 2 ] 1/2 > 



(13) 



where the primed quantities refer to the case without soft- 
ening. 

In general, the value of r max tends to be much more af- 
fected by the softening than Vm ax - This can be clearly seen in 
the upper panel of Fig. [3] where we show the Vm a x-r max rela- 
tion for all main haloes of the MS-II. Note that the relation 
has some intrinsic scatter, but we here only show the mean 
values of r max in logarithmic bins of Vmax- The red solid lines 
show the correlation for different redshifts equal to z = 0.5, 
1.0, 1.5 , 2.1, 3.1, and 4.9, after the softening correction of 
Eqs. ()13[) has been applied, whereas the thicker line is for 
2 = 0. T he black solid line is an extrapolation of the results 
found bv lNeto et ail (120071 ) at z = for the MS-I. There is 
a clear upwards deviation from the power-law behaviour for 
haloes with particle number below N p ~ 3600 (correspond- 
ing to a mass of ~ 2.5 x 10 10 h~ 1 M & , or Vmax ~ 60kms~ 1 
at z = 0). This deviation indicates that the softening cor- 
rection is only able to correct part of the resolution effects. 
We can be certain that this is indeed a numerical effect by 
looking at the corres ponding results for the Aquarius haloes 
(|Springel et alj 2008), which demonstrate with simulations 
of the same initial conditions but drastically different mass 
resolutions that the r max -V max relation is indeed a power 
law over a large dynamic range, at least for Vmax as small as 
~ 1.5km s -1 . This is t he case for both , main haloes and sub- 
haloes (see Fig. 26 of ISpringel et al] l|2008h ). We note that 
haloes with N p ~ 1000 or less are affected importantly by 
two body relaxation. The discrete representation of a dark 
matter halo with this low number of particles will undergo 
two-body encounters whereas a real dark matter system with 
~ 100 GeV particles is essentially collisionless. Generically, 
the two-body relaxatio n process tends to red uce the concen- 
tration of objects (e.g. iDiemand et al.ll2004T ). We think this 
is the main mechanism responsible for the upturn shown in 
Figs. [3] and [7] even though a detailed study of this effect is 
out of the scope of this work. 

Since an overestimate of r max can result in an underesti- 
mate of the computed luminosity for a given halo, we further 
correct the r max values of small haloes phenomenologically 
in order to allow use of all haloes in the MS-II down to the 
minimum resolved halo mass of 1.4 x 10 s /i -1 M0, and to get 
the correct mean halo luminosity over the full mass range. 
For this purpose, we force the values of r max for all haloes 
with less than 3600 particles to follow the power law 



= A(z) 



( V 

I v nit 



(14) 



where A(z) and ct(z) are the parameters of the power law. 
These parameters actually depend very mildly on redshift 
(with A(0) w 0.04 and a(0) « 1.32 at z = 0). The redshift 
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Figure 3. Uppev panel: Relation between r ma x and Vmax for 
all main haloes in the MS-II at different redshifts. The red solid 
lines are based on measurements corrected for softening effects 
with Eqs. J 1 3 D - The blue dashed lines show the values we adopt 
for calculating the expected luminosity of a halo; here an addi- 
tional correction for r max was applied for small haloes such that 
they extend the power law behaviour measured for well resolved 
structure s. The black solid line is an extrapolation of the results 
found bv lNeto etafl [|2007ri for the MS-I, while the black dotted 
line is the power law obeyed by subhaloes, see Fig. The thicker 
lines correspond to z = 0, while the rest of the lines are for red- 
shifts: 0.5, 1.0, 1.5, 2.1, 3.1 and 4.9. Bottom panel: The relation 
between concentration C200 and Vmax. The line styles and colors 
are the same as in the upper panel. 



dependence can be understood as a change in the mass- 
concentration rel ation of haloes, which gets flatte r for in- 
creasing redshift (|Zhao et al.ll2003l : iGao et alj|2008h . We es- 
timate the concent ration of the main h aloes following the 
description given in lSpringel et all l|2008h . where a measure 
of the concentration is obtained in terms of the mean over- 
density within r max relative to the critical density p C rit: 



p(r a 



= 2 



Pcrit 



Vm 



(15) 



H(z)r m: 

where H(z) is the Hubble rate as a function of redshift. 
The overdensity Sv can be related to the usual concentra- 
tion parameter C200 = ?~2oo / rs °f the NFW profile using the 
characteristic overdensity: 

S c = — — ^2° = _Ps_ = 72USv 

3 ln(l + C200) - C/(l + C200) Pcrit 

(16) 

Using Eq. ()16|) we can compute C200 for a given value of 
Vm ax and r max - The lower panel of Fig. [3] shows the relation 
between C200 and Vmax for the corrected and uncorrected 
values of r max as blue dashed and red solid lines, respectively, 
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Figure 4. Total luminosity coming from main haloes per differ- 
ential mass interval as a function of mass. The thin blue lines 
are for different redshifts as in Fig. \3\ while the thick solid line 
shows the z = result. The dashed lines show the run of our 
extrapolations as discussed in the text. 



Figure 5. Flux multiplier f(M^ > M m ; n ) for the main haloes in 
the MS-II as a function of M m i n . The solid-blue and dashed-black 
lines are analogous to the ones in Fig. [4] The solid-black line is a 
theoretical estimate as described in the text. 



with the redshift increasing from top to bottom. The change 
in slope of the power law followed by the blue dashed l ines is 
rou ghly in agreement with the results of IZhao et alj (|2003l ) 



and 



Gao et al.1 (|200ah . 



4.1 Gamma-ray luminosity of haloes down to the 
damping scale limit 

Following the formulation in section 2, we here analyze the 
flux multiplier for large volumes, Eqs. (O and (J5]), for the 
MS-II. Recall, the flux multiplier gives the ratio of the 7- 
ray flux coming from all haloes inside the simulated volume 
with masses larger than a minimum mass M m i n to the emis- 
sion produced by a homogeneous distribution of dark matter 
filling the box of volume Vb with an average density ps- 

We obtain this dimensionless flux multiplier by defining 
first the function 



F h (M h ) = 



M h A log M h ' 



(17) 



where the sum is over all the luminosities Lh of haloes with 
masses in the logarithmic mass range: \ogMh ± AlogM^/2, 
where A log Mh is a fixed logarithmic bin size; Mh is the 
mean value of the halo mass in the given bin. Using this 
definition we can approximate Eq. Q as 



f(M h > A/ min ) 



1 



P%V B 



F h (M h 
In 10 



■ dM h 



(18) 



In this sense, the function Fh(Mh) is just the total lumi- 
nosity of haloes in a mass range, per unit mass range. The 
function Fh(Mh) is shown in Fig. [4] for all main haloes in 
the MS-II. The different blue lines are for different redshifts, 
as in Fig. [3] 



At intermediate mass ranges, Fh(Mh) is clearly well ap- 
proximated by a power law: 

F h {M h ,z)=A h {z)MZ hiz) . (19) 

Our goal is to fit the parameters of this power law so that 
an extrapolation can be done down to the cutoff mass for 
neutralinos. For the neutralino mass corresponding to the 
model we have chose n, the free streaming mass is of the or- 
der of 10" 7 /i _1 M Q (|Hofmann et alj|200ll ). however, acous- 
tic oscillations due to the coupling between cold dark mat- 
ter and the radiation field in the early Universe, can also 
produce a dam ping in the power spectru m of density per- 
turbations (e.g. lLoeb fc Zaldarriagall2005l ). The cutoff mass 
of the smallest haloes that can be formed is determined by 
the stronge s t of t hese effects. Taking the recent results of 
iBringmannl (|2009h (see their Fig. 3), this cutoff mass for 



185 GeV lies in the range 10 



10~ 4 M Q . We will 



take a fiducial value of 10~ B /i~ 1 Mq for our extrapolation, 
noting that the value of the minimum mass for bound neu- 
tralino dark matter haloes is a source of uncertainty in our 
results. 

We obtain the parameters of the power law in Eq. (|19[1 
by fitting the function Fh(Mh) between two mass limits, 
with the lower limit chosen as Mii mjm i n = 6.89 x 10 8 /i _1 Mq, 
corresponding to haloes with 100 particles (below this num- 
ber the mass and abundance of haloes is not reliable), and 
the higher limit set equal to the last logarithmic mass bin 
with more than 500 haloes, such that uncertainties from 
counting statistics are avoided. We find that for these mass 
ranges, the parameters of the power law fits change only 
slightly with redshift; in fact for z < 2.1, ah — —1.05 with 
less than 2% variation, and Ah — 6.92 x 10 with less than 
50% variation. The black dashed lines in Fig. [3] show the re- 
sulting extrapolation of the power law down to 10 7 ft _1 M0. 
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The blue lines in Fig. [S] give results for the flux multi- 
plier, Eq. (| 18 P . as a function of M m i n for the main haloes 
in the MS-II for different redshifts, as in in Fig. [4] The 
black dashed lines show the extrapolation of the flux mul- 
tiplier down to M m in = 1(T 6 h~ 1 M e using Eq. !fl9)l. The 
flu x multiplier has bee n previously computed analytically 
bv lTavlor fc Silk] (|2003T ). We follow their procedure and cal- 
culate an analytical estimate of f(M > M m i n ) for the cos- 
mology corresponding to the MS-II. We first use Eq. ((3]) to 
compute the individual flux multiplier /^(ca) of an NFW 
halo with a concentration ca given by the redshift mass- 
c oncentration relat ion, computed as in the analytical model 
of lEke et all (|200lT ) and with a choice of A = 178 Q° n i5 ~ 95. 
Using the value of /^(ca) we then compute the full flux mul- 
tiplier based on Eq. (El) , assuming a mass f unction calculated 
with the formalism of ISheth et al.l (|200lh ■ The result is in- 
cluded in Fig. [5] with black solid lines. 

We see that these analytic estimates roughly agree with 
the actual numerical results from the MS-II over the resolved 
mass range. However, small differences in the slope of f(M > 
Mmi n ) at the minimum resolved mass range (1.39 — 6.89 x 
10 s h~ 1 M@), produce a large difference for the extrapolated 
values at the damping scale limit, which amounts to a factor 
of up to 3 — 5. We note that for redshifts z > 2.1, our 
extrapolation down to the cutoff mass is not accurate any 
more and the obtained values are clearly overestimated. This 
is a reflection of the difficulties to still reliably fit a power 
law to F(M) at these high redshifts, where the population 
of haloes over the resolved mass range becomes ever smaller. 

Another interesting question that arises from this sta- 
tistical analysis is: Which mass scale contributes most to the 
total flux of gamma rays coming from dark matter annihi- 
lations? The answer is presented in Fig. [5] which shows the 
differential cumulative flux, which is really just the loga- 
rithmic derivative of the cumulative flux shown in Fig. [5l 
as a function of the minimum halo mass. From the fig- 
ure we see that less massive haloes contribute increasingly 
more to the total gamma-ray flux coming from dark mat- 
ter annihilation. There is a clear change in behaviour for 
a given mass scale M m i n ~ 10 14 ft _1 M0 at z = 0, going 
down to Mmin ~ 10 12 h~ 1 M Q at z = 2.1. Fig. [6] also re- 
veals an artificial downturn for the minimum mass range 
resolved (1.39 - 6.89 x 10 8 h _1 M ). The analytical expec- 
tations, shown as black solid lines, qualitatively agree with 
the results from the simulation. However, the quantitative 
differences in the extrapolated region are large. 

4.2 The effects of substructure 

Substructures within CDM haloes produce a total 7-ray 
luminosity which is dominant over t he smooth halo com - 
ponent for an exter nal observer fe.g iTavlor fc S ilk 2003). 
ISpringel et alj (2008) showed that for a MW halo, substruc- 
tures with masses larger than 10 M© within r2oo emit as 
much as 76% of the emission of the smooth halo. Also, using 
their set of simulations of the same halo at different resolu- 
tions they extrapolated the contribution from substructures 
down to masses of 1O _6 M0, at which point they found the 
total cumulative luminosity to be 232 times larger than the 
contribution from the smooth halo component. 

In this subsection, we analyse the contribution of sub- 
structures to the 7-ray flux of the main haloes using as a 
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Figure 6. Differential cumulative flux multiplier |d/(Af n > 
M min )/dlogAf m in I as a function of minimum mass. The line 
styles and colors are as in Fig. [5] 
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Figure 7. Relation between r max and Vmax for the subhaloes of 
the 10 most massive clusters in the MS-II at z = 0. These main 
haloes have Vmax values in the range ~ 850 — 1250kms _1 . Line 
styles and colors are as in Fig. [3] 

basis the most massive haloes in the MS-II that contain 
significant resolved subhaloes and are therefore suitable for 
such an analysis. To this end we follow an analogous proce- 
dure to the one used in the previous subsection. 

Fig. [7]shows the r max — V max relation for the population 
of subhaloes inside the 10 most massive clusters of the MS- 
II at z = 0, the line styles and colors are the same as in 
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Fig. [3] These clusters have Vmax values in the range 850 — 
1250 kms -1 . The subhaloes follow a very similar power-law 
behaviour as the population of main haloes in the upper 
panel of Fig. [3] The slope of the power law is the same, only 
the normalization is slightly different, a factor of 0.62 lower 
than for main haloes. This is in agreement with fin dings 
obtained for the Aquarius project (|Springel et al.ll2008h . 

As for main haloes, we have corrected the values of r max 
and Vmax of the subhaloes for softening effects according to 
Eqs. (|13[) . Also, as for the main haloes, the upturn in the 
fmax — Knax relation for small subhaloes clearly indicates 
an overestimation of r max in this regime due to numerical 
resolution effects. In order to avoid biasing our results, we 
hence force the subhaloes in this low-r- max low- Vmax regime 
to continue to follow the power law defined by the larger, 
more massive subhaloes: 



^ S ub(z)(V m ax) C 



(20) 



where a su b has a slight redshift dependence that we can 
neglect since a value of a su b ~ 1.34 fits the mean behaviour 
of the data for z < 3 with an accuracy better than 20%. We 
find that the value of A su b(z) can be approximated as: 



,(z) ~ 0.023(1 + z) 



-1/2 



(21) 



As a further consistency check, we have analyzed all main 
haloes in the MS-II with more than 500 subhaloes, finding 
that they agree well with the behaviour of the 10 most mas- 
sive clusters shown in Fig. [7] Thus, all resolved subhaloes in 
the MS-II with particle number N p < 3600 were corrected 
using Eq. ((20)) . 

The contribution of subhaloes, per unit mass range, to 
the total 7-ray luminosity of their host can be analyzed using 
the following quantity: 



M sut 
M h 



Mh \ E^sub 

L h ) A/subAlogM sub 



(22) 



where M su b and L BU ^ are the mass and gamma-ray luminos- 
ity of a given subhalo. We have computed F su h for all main 
haloes in the MS-II with more than 500 subhaloes and found 
that this quantity can be described by a power law in the 
intermediate mass range, in a similar way to the function 
Fh for main haloes: 



Fsuh 



~Mh~ 



(23) 



We estimate the parameters of this power law by fitting the 
results between M su b, m m = 6.89 x 10 8 ft _1 M0, the same 
lower mass limit used for main haloes, and an upper mass 
limit given by the most massive logarithmic mass bin with 
more than 15 subhaloes. Using these two mass limits a power 
law is found to provide a good fit to the general behaviour 
in the intermediate mass range for all main halo masses and 
at different redshifts. 

Fig. [5] shows the parameters ^4 S ub and a su b of the fit 
as a function of the mass of the host. The median values 
for logarithmic mass bins at z — are shown with thick 
black solid lines, the first and third quartiles with thiner 
black solid lines. The dark grey and light grey lines are for 
z — 1.0 and z = 2.0, respectively. The results show almost 
no dependence of the parameters of the power law fit on the 
mass of the host and on redshift. For the majority of the 




10' 



10 M 

K Kit'] 

Figure 8. Parameters ct su b and Ajub from Eq. (12.'St as a function 
of main halo mass. The thickest solid lines are the median values 
for logarithmic mass bins, the dotted regions are within the first 
and third quartiles. The black, dark grey and light grey lines are 
for z = 0, 1.0 and 2.0, respectively. 



cases, these parameters lie in the range: 



■0.95 sS a suh < -1.15 
-0.5 «S logyl SU b < 0.1. 



(24) 



For lower host halo masses, the scatter in the parameters 
for a given mass bin becomes increasingly larger, however. 
In these cases, the low number of substructures per halo 
makes the power-law fit less reliable. 

We will assume that the power law in Eq. (I23[) is uni- 
versal, allowing us to extrapolate the value of F su b down to 
the damping scale limit. Using the range in Eqs. (I24p . we 
can obtain estimates for the maximum and minimum con- 
tribution of all substructures to the gamma-ray luminosity 
of their host halo: 



/sub (/max, Mh) 



1 



/maxM,, 
10~ 6 



-1 



In 10 



dM sub 



(25) 

where / ma x is the ratio of the most massive subhalo to the 
mass of the host. If a su b 7^—1 then: 



/sub (/max, Mh) 



(/max) 



my 



In 10(a su b + 1) 

(26) 

For the range of values given in Eqs. (|24[) . which reflects the 
uncertainties in our power law extrapolation, the total lumi- 
nosity of all substructures in a halo of mass 10 12 ft _1 M s lies 
in the range / su b £ (2, 1821) times the luminosity of the main 
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halo, bracketing; the value of 232 found by ISpringel et ail 
(2008) for the Aq-A halo of the Aquarius project. 



5 CONTRIBUTIONS TO THE EGB FROM 
DARK MATTER ANNIHILATION 

With the results of the previous sections, we finally concen- 
trate on the main objective of this work, namely, to calculate 
the contribution of dark matter annihilation to the EGB. 
Specifically, we are interested in making realistic maps of 
the expected specific intensity across the sky. 



5.1 Simulation of the past light cone 

The specific intensity (Eq. 1 1 1 p describes the total emission 
from dark matter annihilations integrated over the full back- 
wards light cone along a certain direction. We use the data 
of the MS-II to fill the whole volume contained in the past- 
light cone of an observer located at a fiducial position in the 
box at z = 0. As we have 68 simulation outputs in total, we 
can approximate the temporal evolution of structure growth 
by using at each redshift along the past light-cone the out- 
put time closest to this epoch. To cover all space, we use 
periodic replication of the simulation box. 

However, this periodic replication would introduce 
strong correlations along the line-of-sight on the scale of 
the box size, and in particular, would lead to replications of 
the same structures along certain sight lines. To avoid this 
problem, we subdivide the backwards light-cone into shells of 
comoving thickness equal to the boxsize. Within each shell, 
an independent random rotation and translation is applied 
to the pattern of boxes that tessellate s the shell, as sketche d 
in Fig. [9] This procedure, first used in lCarbone et al] l|2008h . 
eliminates the unwanted line-of-sight correlations, while at 
the same time it maintains continuity of the maps in the 
transverse direction across the sky. As a result, our final 
maps are free of obvious tessellation artifacts. We note, how- 
ever, that this procedure can not make up for the missing 
large scale power on scales larger than the box size itself. 
Also, it will not be able to eliminate potential excess power 
in the angular power spectrum on the periodicity scale of 
the box, as within each shell the simulation box seamlessly 
tiles the shell. However, the fixed comoving box size trans- 
lates to a different angular scale in each shell, so that the 
averaging over many shells usually makes this a small effect. 

As discussed earlier, our strategy is not to account for 
the emission at the level of individual particles, but instead 
to use entire haloes and subhaloes, allowing us to accu- 
rately correct for resolution effects. The list of dark matter 
particles at each output time is hence replaced by a cata- 
logue of dark matter substructures, each with a known to- 
tal luminosity and known radial luminosity profile, calcu- 
lated from their r max and V max values as described earlier, 
and including a resolution correction for poorly resolved ob- 
jects. Map-making then becomes a task to accurately accu- 
mulate the properly redshifted emission from these struc- 
tures in discretized representations of the sky. For the pix- 
elization of the sky we use equal area pixels based on the 




Figure 9. Sketch of the light-cone reconstruction carried out with 
the periodic MS-II simulation. We subdivide the light-cone into 
shells of comoving thickness 100 /i — 1 Mpc. Within each shell, a 
coherent randomization of the orientation and translation of the 
periodic tessellation with the simulation box is carried out, in 
order to avoid repetition of the same structures along a given 
line-of-sight and at the same time to maintain continuity in the 
transverse direction across the sky. The time evolution of cosmic 
structure is taken into account by using at each position (redshift) 
along a line-of-sight the closest simulation output time. 

HEALPIX0 tessellation (|G6rski et al.l 120051 ). All our maps 
use iVpi x = 12 x 512 2 ~ ir x 10 6 pixels corresponding to an 
angular resolution of ~ 0.115°. 

A given pixel in the simulated maps covers a solid angle 
Af2pi x . Our map-making code computes the average value 
of the specific intensity within the area subtended by this 
solid angle by conservatively distributing the emission of 
each substructure over the appropriate pixels. Combining 
Eqs. ([9]) and (|11|) we can write the relevant sum over all 
substructures in the light-cone as: 

J 7 ,o(Aripix) = — ^2 L h w{d h ,r h )E 1 fifsvSY{z h )\ El0 

heAf2 pix 

(27) 

where the function w(dh,rh) is a weight function that dis- 
tributes the luminosity of a given halo onto the pixels over- 
lapping with the projected "size" (or more precisely the lu- 
minosity profile) of the halo; the latter depends on the dis- 
tance of the observer to the halo du and the transverse dis- 
tance rh between the halo centre and the centre of the pixels 
it touches. Except for structures that are very nearby, the 
high central concentration of the emission of a subhalo and 
the limited angular resolution of our maps give most sub- 
haloes the character of unresolved point sources. It therefore 
makes little difference in practice if the transverse luminos- 
ity profile of a subhalo, a p 2 -weighted projection of a NFW 
profile, is replaced with an SPH-like kernel function in 2D 
with radius equal to the halo's half- mass radius T\/2- How- 
ever, it is important to guarantee that the sum of the weights 

3 The HEALPIX software package is available at 
http://healpix.jpl.nasa.gov 
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w{dh, Th) associated with all pixels that overlap with a given 
substructure 'add up to unity', i.e. that the total emission of 
the halo is represented in the map in a conservative fashion. 
Our map-making code guarantees this property. For perfor- 
mance reasons, the code also uses a tree-based approach to 
quickly find all subhaloes (and also any of their periodic im- 
ages, if appropriate) that overlap with a given line of sight. 

For the purpose of extending the predictions of our 
maps down to the damping scale limit of neutralinos (~ 
1O _6 M ), we will incorporate the results of section 4 by 
dividing the maps into separate components: (i) the con- 
tribution of resolved haloes and subhaloes with a minimum 
mass of 6.89 x 10 s 7l _1 M (called "ReHS" in the following); 
(ii) the contribution of unresolved main haloes with masses 
in the range 1.0 x 10" 6 - 6.89 x 10 s /i _1 M Q (referred to as 
"UnH"); and finally, (iii) the contribution of unresolved sub- 
haloes in the same mass range as in the case (ii) (component 
"UnS"). 



5.2 Resolved haloes and subhaloes (ReHS) 

In Fig. 1 1 01 we show full sky maps of the the 7-ray emission 
of all main haloes and subhaloes that are detected in the 
MS-II; this corresponds to an essentially perfectly complete 
sample of all haloes above a mass limit of 6.89 x 10 s /i _1 M . 
In the top panel of Fig. [T0l a partial map at low redshift (cor- 
responding to the first shell at z = of depth 50/i _1 Mpc) 
is shown in order to illustrate typical foreground structure, 
whereas the bottom panel gives an integrated map out to 
2 = 10, which is approximately the full EGB from annihi- 
lation. The maps in Fig. [TU] are for _E 7i o = 10 GeV, where 
E-yfi is the photon energy measured by the observer (i.e. the 
emission energy of the photons is _E 7 = (1 + z)E 1 ,o). In the 
color scale used for the maps, red corresponds to the high- 
est and black to the lowest values of the specific intensity, 
respectively. 



5.3 Unresolved haloes (UnH) 

The maps shown in Fig. \W\ are complete down to the 
minimum mass in the MS-II that we can trust, ~ 6.89 x 
10 8 /i _1 M . In order to make a prediction of the full EGB 
coming from dark matter annihilations, we extrapolate the 
7-ray flux to account for the contribution of all missing dark 
matter haloes down to the cutoff mass 10~ 6 fc _1 M s . For 
this purpose we use the results found in section 4 on the 
universal power law behaviour of the functions Fh(Mh) and 

-fsub 

(M Buh /M h ). 

The way we incorporate this contribution in the 7-ray 
maps is the following. We assume that the EGB radiation 
from the missing haloes in the mass range 10 -6 /i _1 M to 
~ 6.89 x 10 8 /i _1 M is distributed on the sky in the same 
way as the one from the smallest masses we can resolve in 
the simulation, which we adopt as the mass range between 

1.4 x 10 8 /i _1 M Q and ~ 6.89 x 10 8 /i _1 M Q (haloes with 20 
to 100 particles). This assumption can be justified by the 
clustering bias of dark matter haloes that appears to ap- 
proach an asymptotic const ant value for low halo masses (e.g 
iBovlan-Kolchin et alj|200gh . Hence, we compute the value 
of a boost factor bh with which each resolved halo in the 
mass range 1.4 — 6.89 x 10 s /i _1 M needs to be multiplied 



such that the luminosity of the unresolved main haloes is 
accounted for as well. Using the power law behaviour of 
F h {M h ) and f(M h > M h , min ), Eqs. (QU and JT5J, we find 
that: 



b h = 



/(kt 6 hr^u^m x 10 8 ft^Mp). 

/(1.4 x 10 8 ft- 1 M ,6.89 x 10 s h^Ms)^, 



60, 



(28) 

where bh is given by the ratio of the flux multipliers com- 
puted analytically, between the cutoff mass limit and the 
100 particle limit, and computed within the simulation in 
the last resolved mass range (20-100 particles). We find that 
bh is nearly independent of redshift up to the highest red- 
shift we can reliable measure the power law Fh(Mh), that is 
for 2 < 2.1. 



5.4 Unresolved subhaloes (UnS) 

To add the contribution of unresolved substructures to the 
7-ray maps, we first assume that the identified subhalo 
population is complete down to a mass limit of M su b = 
6.89 x 10 8 ft _1 M . Below this mass we compute luminosi- 
ties using the formulas presented in section 4.2, separating 
the way we add this contribution in two cases. 

In the case of resolved main haloes with more than 100 
particles, we simply use Eq. (|26[) with / max Mh = 6.89 x 
10 8 /i _1 M if the main halo has subhaloes, and / max = 0.05 
otherwise. For the former case we distribute the extra lu- 
minosity coming from unresolved substructures among the 
subhaloes of the main host, assuming in this way that un- 
resolved subhaloes are distributed in the same way as the 
resolved ones. For the latter, the extra luminosity is given to 
the main background subhalo. We note that in this case the 
particular value of / max has little impact in Eq. (126f) since 
M h > 6.89 x 10 8 ft _1 M and thus, the right-hand side of 
the term in brackets is dominant for a su b < — 1; for the ex- 
treme value of a su b = —0.95, values of / max between 0.005 
and 0.5 give variations of / su b smaller than 3.5%. We use 
/max = 0.05 since this is a typical value in the MS-II. 

For the second case of all main haloes, resolved and un- 
resolved, with masses less than 6.89 x 10 s /i _1 M we use 
the following strategy. Since in the last subsection we de- 
scribed how to get the total luminosity coming from all 
main haloes with masses between the damping scale limit 
and 6.89 x 10 8 /i _1 M , we can compute the boost factor 
fe su b for which this luminosity needs to be multiplied in or- 
der to include the full contribution of all subhaloes of the 
main haloes in this mass range: 



k 

where 



fboos t (lQ- 6 ft"'M s , 6.89 x 10 8 /t^Mp) 
mb fno-boost(10- 6 fc-iMg, 6.89 x 10« fc-iM ) [ ' 



, : no-boost is given by: 
jno-boo St (10- 6 /i _1 M ,6.89 x 10 s ft _1 M Q ) 



Fh(M h 



Jio-6 In 10 

and fboost can be written as: 
fboo S t(10 _6 /i _1 M Q ,6.89 x 10 s r'Me) 



dM h 



(30) 



/•6.89x10 s 



>10~ e 



, f ( f M y Fh(Mh) 

+ hub (,,/max, Mh )\ — j^j — 
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Figure 10. Upper panel: One of the partial maps (z = 0) showing the cosmic 7-ray background produced by dark matter annihilation. 
The color scale gives a visual impression of the values of the specific intensity for each pixel in the map; the red color corresponds to 
the highest values of specific intensity. The observed energy of the simulated 7-ray radiation is lOGeV, and the benchmark point L as 
described on Table 1 was used as input for the supersymmetric model. Lower panel: Co-added map showing the full 7-ray sky map from 
dark matter annihilation integrated out to 2 = 10. 
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After using Eqs. (|18[1 and (|19[1 we can write: 

6 sub = i + . m( Au \ J f2s? +1 (io- 6 )— +1 

lnlO(a sub + 1) I 

\a h - a sub / (6.89 x lO 8 )"^ 1 - (lO" 6 )"'^ 1 J ^ ^ 

We found that 6 su b is roughly independent of redshift. Using 
the values of ah = —1.05, A B ub, a su b determined earlier, we 
obtain 6 su b ~ 2 — 60. We shall take the extreme values of this 
range in the following presentation of our results, as these 
minimum and maximum values reflect the uncertainties in 
our extrapolation and should bracket the true result. 



5.5 Isotropic and anisotropic components of the 
EGB from dark matter annihilation 

The upper panel of Fig. [11] shows the mean value of the 
specific intensity A/ 7 .o divided by the comoving thickness 
of the shell (equal to 100 ft, _1 Mpc, except for the map at 
2 = where the shell is half as thick) at i? T ,o — 10 GeV for 
our partial sky maps as a function of redshift. Recall that 
each partial map represents the total specific intensity in a 
shell of constant comoving thickness, so this provides infor- 
mation about where the signal is coming from. The black 
line is for the ReHS case, while the dashed line is for main 
haloes only (resolved and unresolved, including the boost 
factor bh = 60). The remaining lines are for the resolved 
and unresolved main haloes and their subhaloes, boosted 
to include the contribution of structures all the way down 
to the damping scale limit; here the light grey line is for 
the "minimum boost" with b su b = 2 whereas the dark grey 
line is for the "maximum boost" with 6 sub = 60, but both 
include the boost bh = 60 for unresolved main haloes. 

In all cases, the contribution from partial maps up un- 
til 2 ~ 1 is almost constant, indicating that the increasing 
number of sources for shells at higher redshifts, due to the 
larger volume seen behind each pixel, is approximately com- 
pensated by the distance factor and the spectral effects from 
the cosmological redshifting and the intrinsic variation of the 
emission spectrum, as described by /susy (see Fig. [32J. 

The lower panel of Fig. [TT] shows the accumulated mean 
value of the specific intensity as a function of redshift. This 
clearly shows that the most relevant contributions come 
from maps up to z ~ 2. For higher redshifts, the partial 
maps contribute less than 5% to the total. Overall, applying 
the maximum boost increases the mean value of J 7 ,o of the 
ReHS case by three orders of magnitude. However, it is im- 
portant to note that Fig. Illl depends on the observed energy 
i? 7 ,o, and more specifically on the spectral shape of the fac- 
tor /susy in the energy range E 1}0 - -£7,0(1 + 2 max ), where 
2 m ax is the largest redshift considered, in our case 2 max = 10. 
For most energies, the /susy spectrum is monotonically de- 
creasing. In these regions, the dependence on 7 7i o with red- 
shift will have the generic form shown in Fig. 1111 However, 
at the highest energies the shape is expected to change dra- 
matically due to the importance of Internal Bremsstrahlung 
and/or of monochromatic lines (see Fig. [2} close to the rest 
mass energy of the dark matter particle. Such an effect is 
clearly visible in the black dotted line in Fig. 111! which is the 
ReHS case for an energy of 100 GeV and features a bump 




1+z 

Figure 11. Upper panel: Mean annihilation intensity per comov- 
ing shell thickness, AJ 7 o/Ax, as a function of redshift for indi- 
vidual partial maps. The solid black line is for resolved haloes 
and subhaloes. The dashed black line is for main haloes only, 
resolved and unresolved down to the cutoff mass using a boost 
factor = 60. The light grey and dark grey solid lines show the 
contribution of all components, resolved and unresolved haloes 
and subhaloes, boosted with the extreme values in the inter- 
val fe su b = 2 — 60. All lines are for -E 7 ,o = 10 GeV, except for 
the black dotted line which is for the ReHS component with 
^7,0 = 100 GeV. Lower panel: The same as the upper panel 
but for the accumulated intensity, 7 7l o(> 1 + z). 



in the redshift distribution, a reflection of the importance of 
IB at the highest energies. 

The angular power spectrum of the EGB is an impor- 
tant tool to study the statistical properties of its anisotropy 
on the sky. In fact, certain classes of 7-ray sources are ex- 
pected to exhibit different power spectra, making this a po- 
tential means to identify the origin of the unresolved EGB. 
If we decompose the relative deviations of the EGB around 
its mean into spherical harmonics with coefficients a; m , 

a^m) = j7 ' o(g ;f ) ~ (/7 '" ) = £ 2 ««.*wm>. 

[1 ^' 0> l = 0m=-l 

(33) 

then the angular power spectrum is defined by the real num- 
bers 

C ' = 2m(i:W 2 )' ( 34 ) 

For small angular scales (i.e. large I), the variance of Aj 
per unit decade in I can be connected directly to the C;'s 
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Figure 12. Angular power spectrum of the EGB produced by 
dark matter annihilation as a function of the multipole I for 
E-yfi = 10 GeV. The solid black line is for resolved haloes and sub- 
haloes. The dashed black line is for main haloes only, resolved and 
unresolved down to the cutoff mass using a boost factor = 60. 
The light grey and dark grey solid lines show the contribution 
of all components, resolved and unresolved haloes and subhaloes, 
boosted with the extreme values in the int erval b sub = 2 — 60 . 
The dotted lines show the predictions from lAndo et af. | <2007bf l 
for a subhalo- dominant contribution with and without consider- 
ing tidal destruction (lower and upper dotted lines), respectively. 
The dash-dotted line shows the case where the main halo domi- 
nates the signal according to the same study. 



through 



d < A U ,,,y 



(35) 



In Fig. 1121 we show the angular power spectrum 1(1 + 
l)Ci/2ir as a function of the multipole I, at an observed en- 
ergy of E 1 fi — 10 GeV, for all the cases discussed in Fig. 1111 
At large scales, I < 10, the power spectrum is related to the 
clustering of dark matter haloes. All cases have the same 
shape on these scales. However, when only the main halo 
contribution to the EGB is considered (dashed black line), 
the normalization is lower because most of the 7-ray signal 
comes from low mass haloes that are less clustered (biased) 
than more massive haloes; recall that in this case we extrap- 
olate the signal down to the cutoff mass by assuming that 
the unresolved main haloes are clustered in the same way as 
the least massive resolved haloes in the MS-II. 

The light grey line in Fig. 1121 corresponding to the 
full extrapolation including subhaloes but with the minimal 
boost & su b = 2, has exactly the same power spectrum than 
the case with main haloes only, at all scales. This is because 
the signal from main haloes is dominant in this case, and 
subhaloes have a negligible effect. In contrast, in the case 
where subhaloes have a significant contribution mediated 
by 6 su b = 60 (dark grey line), the normalization is larger at 
small scales because the signal is dominated by subhaloes 



belonging to the most massive haloes, and the latter are 
strongly clustered. 

As the angular scale decreases, I > 10, the power spec- 
trum depends more and more on the internal structures of 
haloes. For the cases where substructures are ignored or are 
negligible (dashed-black and light grey lines), the slope be- 
comes steeper, with a slope close to 2. The power spectrum 
for the cases where substructures are relevant (dark grey and 
black solid lines) behaves differently, however. In the range 
I £ [20, 100], it becomes slightly shallower, i.e. the signal is 
slightly more isotropic in this regime. This is probably pro- 
duced by the distribution of substructure within the biggest 
haloes. Contrary to the strong central concentration of the 
matter in a halo, the number density profile of subhaloes 
is considerably shallower than a NFW profile and produces 
a luminosity profile in projection which is essentially flat 
|Springel et ail 120081 ). This effect continues until I ~ 200 
where the power spectrum becomes dominated by the low- 
mass main haloes. The corresponding upturn happens at 
smaller scales for the case corresponding to the black line 
because the substructures in the most massive haloes re- 
main important for a larger angular range, whereas for the 
case of the dark grey line, the low-mass unresolved main 
haloes start to dominate earlier. 

The dash-dotted and dotted l i nes in Fig. [T2]were taken 
from the analysis of lAndo et ail l|2007bl ). They were ob- 
taine d under the assum ption of a Sheth-Tormen mass func- 
ti on (ISheth et al.ll200lh . a halo bias according to the model 
of lMo fc White! (|l996T ), and a NFW main halo density profile 
with a conce ntration given by the mass-concentration rela- 
tion found bv lBullock et ail l|200ll) . Subhaloes were incorpo- 
rated through a halo occupation distribution model, where 
the number of subhaloes in each main halo, (N su h\Mh) , is 
assumed to depend only on the mass of the host according 
to (N sub \M h ) oc Ml. The upper dotted line in Fig. 1121 is 
for P — 1.0 and the lower dotted line for f3 = 0.7; in these 
cases subhaloes dominate completely the signal. The sub- 
haloes are assumed to be distributed with a NFW radial 
profile inside the main halo. Since they dominate the sig- 
nal, the total gamma-ray luminosity (halo + subhaloes) is 
effectively tracing the density of dark matter instead of the 
density squared as expected for a smooth halo with no sub- 
structures. This case is represented with a dash-dotted line 
in Fig. 1121 In all three cases only the contribution of haloes 
down to a minimum mass of 10 6 Mq is taken in to account. 

Th ere are several approximations used by lAndo et ail 
|2007bT ) that we believe may be the cause of the differences 
between their results and ours. In our case, the clustering 
of dark matter haloes is given directly by the simulations 
instead of by an approximate halo model. This introduces 
differences that can only be quantified by a dir ect compar- 
ison of the power spectrum of haloes used by lAndo et ail 
(|2007br ) and the one in the MS-II. Such a comparison is out 
of the scope of this paper. We need to mention however, that 
our simulated maps have a deficiency of power at the largest 
angular scales due to the lack of comoving scales larger than 
the simulation box (100 Mpc/h ) . This is a problem already 
pointed out by ICarbone et ail (|2008l ) that used a similar 
map- making procedure, and is probabl y the main reason for 
the disagreement between the results of lAndo et ail l|2007bl ) 
and ours on the large scale regime. 

On the other hand, while we do assume a NFW pro- 
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file, the concentration of each halo is given by its scal- 
ing pro perties (see Eqs. [15] and I16[) , whereas lAndo et al.l 
(2007b) uses the mass-c oncentration fitting relations found 
by iBullock et all (|200lh . This fitting formulae agree well 
with simulations at z = but strongly disagree at higher 
redshift. They predict a relatively strong correlation between 
concentration and mass, contrary to simulations where this 
correlation gets flatter for higher redshifts. This difference 
co ntributes to t he disc repancy between our results and those 
of lAndo et all (|2007bh in the small scale regime. Finally, an 
important disagreement is ex pected due to the tre atment of 
subhaloes. The model used bv lAndo et all (|2007bl ) produces 
a shallower power spectrum at small scales since the gamma- 
ray luminosity of a halo is effectively given by the density 
of dark matter instead of by the square of the density. In 
simulations, the radial distribution of subhaloes is consider- 
ably shallower toward s the centre than a NFW profile due 
to ti dal striping (e.g. iDiemand et all 120081 ; ISpringel et all 
2008). Near the halo centre our expectation is that the sig- 
nal will be dominated by the smooth component rather than 
by substructures. This is an effect that contributes to make 
our power spect rum steeper at the smaller scales than the 
one obtained by lAndo et all l|2007bl ). 

5.6 Energy dependence 

So far we have restricted our analysis to an observed 7- 
ray energy of 10 GeV. In the following we study how the 
isotropic and anisotropic components of the 7-ray maps 
change with photon energy. The black solid line in Fig. [13] 
shows the energy spectrum of the mean annihilation in- 
tensity I 7i o for the case of the maximum boost bu = 60, 
fesub = 60 (corresponding to the dark grey line in Fig. Ill[) . 
For simplicity we refrain from showing the other cases explic- 
itly as they differ only in their normalization (for example, 
the minimum boost case with b su b = 2 simply lies 2 orders 
of magnitude lower) . The solid circles with error bars are ob- 
servational measurements taken by the satellite EGRET a f- 
ter removal of the galactic foregrounds (|Strong et al.ll200^ ). 
There are several sources other than DM annihilation that 
contribute to the EGB in the range of energies shown in 
Fig. 1131 blazars and cosmological structure formation shocks 
are among the most important ones. In Fig. 1131 we include 
the contri bution to the EGB from the blazar model intro- 
duced by iKneiske fe Mannheim! j2008) (dashed line) that 
explains most of the signal for energies lower than 1 GeV. 
However, there seems to be an excess of 7-ray flux observed 
in the EGB for energies larger than 1 GeV that ca n not 
be re adily explained by known astrophysical sources l|Onel 
2006). Interestingly, the 7-ray spectrum produced by dark 
matter annihilation has the right shape to fit the observed 
values in this energy range, even though this could of course 
be a coincidence. Our maximum boost scenario is in any 
case still too low (by about an order of magnitude) to ex- 
plain the dominant component of the EGB in this energy 
range, and an additional boost factor would be needed to 
explain this signal, if confirmed, with DM annihilation. An 
improved and more precise measurement of the energy spec- 
trum will be made by FERMI very soon, so there is hope 
that this issue will become clearer in the near future. We 
note that completing and improving the measurements for 
high energies -E 7 ,o > 20 GeV is of key importance since the 
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Figure 13. Mean annihilation intensity -E 7 , 0^7,0 as a function 
of the observed energy for the maximum boost case, = 60, 
£>sub = 60. The observational values an d the error bars com e from 
the EGRET satellite, as presented bvlStrong et alj Efloj ). The 
"blazar model" given by IKneiske fc Mannh eim ( 200f|) is shown 
as a dashed line. 

DM annihilation spectrum is expected to drop off sharply at 
such higher energies, making this an important discriminant 
against other types of sources. 

The presence of a secondary bump at the highest ener- 
gies, -E 7 ,o ~ 125 GeV is a consequence of IB effects that be- 
come important at these energies. The relevance of this effect 
depends on the specific parameters chosen for the mSUGRA 
model (see Fig. [5J, but its detection in the isotropic com- 
ponent of the EGB is generically less likely than identifying 
the primary bump. 

Fig. [14] shows the ratio of the power spectrum Ci for 
different energies to the power spectrum at 10 GeV as a 
function of the multipole The lines are labeled in the left 
of the figure according to the energy of the map in GeV. 
For E-yfi £ (0.1, 30) GeV there is generally an increase of 
po wer for increasing energy , a feature already pointed out 
bv lAndo fc Komatsul (l2006t ). Whereas overall the variations 
of the ratio Cj/Cio for these low-intermediate energies are 
small, there are some differences at intermediate to small 
scales that have important consequences, as we discuss be- 
low. For the higher energies, -E 7 ,o £ (100, 156) GeV, the sig- 
nal originates completely at relatively low redshifts. For ex- 
ample at E-yfi = 100 GeV, all the signal comes from z < 1.0 
(see Fig. Hip since the /susy spectrum has a sharp cut-off 
at the neutralino mass, 185 GeV. The contribution of IB to 
the photon spectrum is dominant at these energies. 

We note that there is a prominent peak in the power 
spectra for energies 100 and 125 GeV, at / = 100 and I = 60, 
respectively. This peak is in fact a spurious effect related to 
the periodicity across the sky that is present within a given 
shell (partial map) . Since we are working with only one sim- 
ulation (MS-II) with a finite box size of 100 /i _1 Mpc, we are 
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Figure 14. Ratio of the angular power spectrum of the EGB 
produced by dark matter annihilation at different energies to the 
one at E 7 ,o = lOGeV (Fig. [T2j . The lines are labeled according 
to their respective energy in GeV. They appear in a grey scale 
going from black to light grey, and with decreasing thickness, 
for increasing energy. The line styles alternate between solid and 
dashed for consecutive energies. 

inevitably affected by the periodicity of the replicas within 
a given shell. For a particular shell, this periodicity intro- 
duces a mild correlation at an angular scale determined by 
the comoving distance from the observer to the shell and 
the transverse box size of 100 /t -1 Mpc. This correlation is 
largely averaged out for the full line of sight integration, as 
the relevant angular scale changes with redshift. However, 
for the highest energies, the bump in the spectrum of /susy 
maximizes the contribution from a particular emitted en- 
ergy, which in turn translates into a particular redshift (see 
Fig. that contributes most prominently. As a result, the 
angular correlation from the periodicity corresponding to 
this particular redshift is particularly strong, and remains 
visible in the final power spectrum. 

To make this effect more evident we show in Fig. I15l the 
angular power spectrum of different partial maps (i.e. single 
shells of comoving thickness 100/i~ 1 Mpc along the line of 
sight), for an observed energy of _E 7i o = 10 GeV (to show 
the periodicity effect, the choice of energy is not important). 
The black solid, dark grey dotted, medium grey dashed and 
light grey dash-dotted lines are for shells at z = 0, 0.5, 
1.0, and 1.5, respectively. Clearly, there are peaks at certain 
values of I for the red, green, and blue lines. The values of 
the multipole I where the peaks occur correspond to angular 
scales that are closely related to the boxsize at the distance 
of the shells. Specifically we find 104.4, 116.7, 99.9 /i _1 Mpc 
for 0.5, 1.0, and 1.5, respectively. This effect is not present 
at z = 0, since the first shell, with a depth of 50 ft~ 1 Mpc, 
encloses a volume too small to contain the full box. 

Another important result that we obtain from the anal- 
ysis of the power spectra is that there are differences in 




10" 5 t ■-<" -i i i 

1 10 100 1000 
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Figure 15. Angular power spectrum of the EGB for different 
partial maps, i.e. for the volume contained in a shell of comov- 
ing thickness lOO/i^Mpc at a given redshift, for the sky map 
at observed energy -E 7 ,o = 10 GeV. The black solid, dark grey 
dotted, medium grey dashed and light grey dash-dotted lines are 
for 2 = 0, 0.5, 1.0, and 1.5 respectively. 

the shape of the power spectra in Fig. [14] for different en- 
ergies. We can exploit these differences and enhance the 
power at small or large angular scales by taking the ratio 
of sky maps at different energies. In particular, the case of 
the energies -E 7l o = 0.1 GeV and E-,.o = 32 GeV (black-solid 
and red-dotted lines, respectively) is notable since these en- 
ergies are free from periodicity effects and when combined 
they can enhance the intermediate to large angular scales. 
The upper panel of Fig. [T5] shows the total sky maps for 
Ejfi — 0.1 GeV and E lt o — 32 GeV in the left and right, re- 
spectively, smoothed with a Gaussian beam with a FWHM 
of 5°. Both maps show only small anisotropies, with the ex- 
ception of a couple of prominent structures. The sky map in 
the lower left panel is the ratio of the two maps; it clearly 
enhances the signal of nearby structures. This is even more 
convincing by looking at the nearest partial map at z = for 
E-y ,o = 0.1 GeV shown in the lower right corner of Fig. 1161 
Most of the prominent dark matter structures that can be 
seen in this map are also clearly present in the map with 
the energy ratio. We conclude that the 7-ray sky maps that 
the FERMI satellite will obtain at different energies could 
be used to construct difference maps which may then show 
enhanced correlations with nearby cosmic large-scale struc- 
ture. 

In our results, we have neglected the interactions 
between the gamma-ray photons produced by annihi- 
lation and the lower energy photons from intervening 
starlight produced in galaxies. This interaction enters ef- 
fectively as an absorption term in the specific inten- 
sity in Eq. (fTTT) due to pair production, and is usually 
parameterised by an e xponential factor e ~ t ( z -- e t-,q) (e.g. 
lAndo fc Komatsull2006l ; IFornasa et alll2009r i . iFornasa et all 
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Figure 16. Upper panel: Total sky maps at energies E 7 ,o = 0.1 GeV and E 7l o = 32GeV in the left and right, respectively. The maps 
were smoothed with a Gaussian beam with a FWHM of 5°. Lower panel: Ratio of the maps in the upper panel (left) and the partial 
map at z = for £ 7 ,o = 0.1 GeV (right). 



(2009) give an appr oximate expression for the opacity r fol- 
lowing the work of ISalamon Sz Steckerl (|l998h : t(z,B 7j o) = 
z/3.3(-E 7 ,o/10GeV) -0 - 8 . The absorption term is typically 
small for energies < lOOGeV. We carried out the exercise 
of producing a pair of sky maps for energies of 10 and 100 
GeV incl uding the absorption term following the parameter- 
isation of lFornasa et alj l|2009l ). For E 1 ,o = lOGeV we found 
no significant changes compared to the map without absorp- 
tion; the total value of 7 7> o drops by less than 10% and the 
normalization of the angular power spectrum increases by 
~ 20% preserving the same shape. For E Jt o = lOOGeV the 
differences are more relevant, the mean specific intensity de- 
creases by 50%, the bump shown in the top panel of Fig. 1111 
(dotted line) becomes less pronounced and therefore the an- 
gular power spectrum gets closer in shape to the ones at 
lower energies; this is an expected effect since absorption is 
stronger at higher redshifts and thus reduces the contribu- 
tion of the "energy bump" at this particular map. 

5.7 Additional boost factor 

The role and importance of annihilations of neutralinos for 
the EGB is subject to many uncertainties, related to both 
the intrinsic nature of neutralinos (as a supersymmetric par- 
ticle) and to its spatial distribution (as dark matter) in the 
different structures in the Universe. The former uncertain- 
ties are ultimately connected to the assumed supersymmet- 
ric model, and within the mSUGRA paradigm, are approx- 



imately constrained by the results shown in Fig. [2] In our 
calculations, they are encapsulated in the value of the quan- 
tity /susy- We chose a particular point, 'L', in the allowed 
mSUGRA parameter space as a fiducial case to obtain quan- 
titative results. This benchmark point is close to be the up- 
per bound for the value of /susy and hence represents a best 
case scenario for a "high" value of the annihilation signal. 
This is important in connection with the results obtained for 
the expected specific intensity spectrum (Fig.ll3[). where the 
value for each energy is given by the isotropic component of 
our simulated sky maps. We mentioned earlier that in order 
to match the observational results coming from the EGRET 
satellite in the 1 — 10 GeV region, a further "boost" fac- 
tor of ~ 10 appears necessary. Since the SUSY benchmark 
point L is close to an upper limit to the value of /susy, this 
boost factor can probably not be accounted for by a differ- 
ent choice of parameters in the mSUGRA parameter space. 
Of course, the situation could be different in broader classes 
of SUSY theories outside the mSUGRA paradigm, but is 
certainly true for our analysis. 

We note that there are a number effects we have not ac- 
counted for but which could possibly give such an additional 
boost factor: (i) We have not incorporated the effects of 
forming luminous galaxies onto the dark haloes; it is gener- 
ally believed that the formation of a baryonic disk contracts 
the halo, pulling matter inwards t owards the centre (e.g. 
iMo et al.lfl99sl : iGnedin et alll2004l ). yielding a higher con- 
centration of the dark matter and a boost of the annihilation 
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flux; (ii) supermassive black holes co uld produce "spikes" of 
dark matter in the gala ctic centre (jGondolo fc Silkl 1 19991 ; 
iBertone fc Merritd fcOQS). Such spikes may also yield an en- 
hanced annihilation signal, (iii) mini-spikes of dark matter 
are expected to form around intermediate-mass black holes 
(IMBHs) (with masses in the range 10 2 < M/M Q < 10 6 ); in 
certain scenarios, the gamma-ray luminosity of each mini- 
spike could be as large as the luminosity of the entire halo. 
For a typical Milky Way halo, the numb er of IMBHs is ex- 
pected to be ~ 100 (|Bertone et al.l 120051 ); ( iv) the non lin- 
ear collapse of dark haloes for cold collisionl ess dark mat- 
ter leads to the formation of caustics (e.g. iHoganl l200ll ; 
I White fc Vogelsber g er|[2009l ). which, due to their high den- 
sity can lead to a la rge annihilation ra t e and gamma-ray 
luminosity. However. IVogelsberger et all l|2009l ). using sim- 
ulations of isolated haloes under self-similar and spherically 
symmetric initial conditions, have recently found that the 
contribution of caustics to the annihilation flux is less than 
previously expected. Based on their results, caustics are ex- 
pected to produce at most a global annihilation flux com- 
parable to the one from the smooth halo component; (v) 
finally, the so-called Sommerfeld enhancement, a quantum 
mechanical focusing effect, may lead to an effective increase 
of the cross section for dark matter with low velocity dis- 
persion, such an enhancement could lead to a boost factor 
of sev eral orders of magnitude (jLattanzi fc Silkll2009l ; iBovvl 
I2009D . 



6 CONCLUSIONS 

Neutralinos in supersymmetric extensions of the standard 
model of particle physics are one of the most promising can- 
didates for being the dominant component of dark matter in 
the Universe. One attractive method for their detection is to 
look for the products of their mutual annihilations. Among 
such products are 7-rays that could be produced abundantly 
enough in the dense centres of dark matter haloes to be sus- 
ceptible of detection. In this work, we have studied the con- 
tribution of these 7-rays to the EGB radiation. Excluding 
the contributions from the Milky Way halo itself, all dark 
matter haloes and their substructures within the past light 
cone of an observer in our galaxy contribute to the EGB. The 
main goal of this paper has been to calculate the expected 
emission from this light-cone on the basis of the state-of-the- 
art Millennium-II simulatio n of cosmic structure formation 
|Bovlan-Kolchin et al.ll2009T ). 

The 7-ray luminosity of each halo (or subhalo) in the 
simulation can be conveniently separated into two factors: 
(i) the SUSY factor /susy, encompassing the intrinsic prop- 
erties of neutralinos as SUSY particles (like its mass, cross 
section and photon yield), and (ii) the astrophysical factor 
related to the density squared dependence of the annihila- 
tion rate and therefore to the clustering of neutralinos as 
dark matter. 

We have discussed in some detail the first of these fac- 
tors in Section 2, with the aim of highlighting the energy de- 
pendence of /susy for a particular set of benchmark points 
within the mSUGRA framework. We also discussed the gen- 
eral mechanisms of photon production: decay of neutral pi- 
ons produced during the hadronization of the primary anni- 
hilation products; (ii) loop-suppressed two-body final states 



containing photons; and (iii) internal bremsstrahlung (IB). 
Among the different benchmark points analyzed we found an 
uncertainty of three orders of magnitude in the normaliza- 
tion of /susy- We chose, for definiteness, one of the bench- 
mark points with a high photon yield and a neutralino mass 
of m x = 185 GeV for our subsequent analysis. 

We then analyzed the astrophysical factor in terms of 
the signals expected from main haloes and their embedded 
subhaloes. In order to avoid strong resolution effects, the in- 
ternal structure of haloes was approximated by spherically 
symmetric NFW density profiles, allowing an analytic in- 
tegration of the total luminosity expected from each halo. 
This effectively yields a scaling law where the flux from a 
halo can be determined by two of its structural properties 
alone, U max and r max . Using this scaling law together with 
the abundance of haloes per unit mass range as given by 
the MS-II simulation, we obtained a power law description 
for the contribution of main haloes of a given mass inter- 
val to the total 7-ray luminosity in a cosmic volume. An 
extrapolation of this power law to the damping scale limit 
of neutralinos in the ~ 100 GeV mass range allowed us to 
include the contribution of unresolved main haloes to the 
total gamma-ray luminosity. A similar universal power law 
was established for describing the contribution of subhaloes 
to the luminosity of their host by analyzing the most massive 
haloes in the MS-II containing at least 500 subhaloes. Using 
reasonably conservative error bounds for the large extrap- 
olation undertaken here, we found that the "boost" factor 
due to substructures lies somewhere in the range 2 to ~ 1800 
for a MW-size halo. With these results in mind, we note that 
the uncertainty in the absolute normalization of the astro- 
physical factor in the annihilation contribution to the EGB 
amounts to two to three orders of magnitude. 

We implemented a new map-making procedure that for 
the first time allowed us to construct realistic realizations of 
the expected extragalactic 7-ray sky from dark matter anni- 
hilation. The angular resolution of the maps was chosen to 
be close to the angular resolution of the recently launched 
FERMI satellite, ~ 0.115°. We found that structures up to 
z ~ 2 contribute significantly in the low to intermediate 
range of the energy spectrum. In particular, for an observed 
7-ray energy of 10 GeV, dark matter structures beyond z ~ 2 
contribute less than 5% to the total flux. At the highest en- 
ergies, between 100 — 185 GeV for the particular benchmark 
point we chose, the IB contribution to the photon yield is 
dominant, resulting in EGB contributions that are predom- 
inantly coming from a narrow redshift range, followed by a 
sharp cut-off towards higher redshifts that reflects the cut- 
off in the energy spectrum of the /susy factor. 

Adopting maximum values for the expected signal boost 
due to substructures, we found that our prediction for the 
energy spectrum of the isotropic component of the EGB 
from dark matter annihilations lies approximately one order 
of magnitude below the observed values of the EGB accord- 
ing to the satellite EGRET, after foreground removal, in 
the energy range 1 — 20 GeV, where an apparent excess of 
7-rays has been reported. The latter has been interpreted 
by several authors in the past as a possible signal of dark 
matter annihilation. At the highest energies, IB produces 
a small bump in the energy spectrum that could in prin- 
ciple be observable but that lies orders of magnitude lower 
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than the expected background radiation for the particular 
benchmark point analyzed. 

We studied the anisotropic component of the EGB by 
computing the angular power spectrum of the simulated 
maps. This yielded specific predictions for the shape of the 
power spectrum, which can potentially be used to discrim- 
inate against other sources of 7-rays, since the annihilation 
signal depends in a specific and unique way on the large 
scale distribution of haloes, on the distribution of subhaloes 
within haloes, and on the abundance and internal structure 
of haloes as a function of time. The shape of the power spec- 
trum was found to depend on the energy of the observations. 
Interestingly, these differences can be exploited to construct 
'color' maps that enhance the signal of nearby dark matter 
structures, akin to hardness ratio maps in X-ray observa- 
tions. For example, we found that taking the ratio of the 
maps at energies of 0.1 GeV and 32 GeV greatly enhances 
the contrast of local dark matter structures, making them 
clearly visible in the 7-ray sky. If strong spectral features 
in the rest-frame emission spectrum of the annihilation ra- 
diation are present, this could be especially powerful, per- 
haps even allowing tomographic observations of dark matter 
structures. 

Although the subtraction of foreground radiation and 
other astrophysical sources contributing to the EGB is a dif- 
ficult task, the specific shape of the contribution of dark mat- 
ter annihilation to the power spectrum of the EGB makes 
this statistical tool attractive for potentially detecting dark 
matter i n the near future (as has been pointed out previously 
by, e.g-. lAndo fc Komatsu|[2006l V The maps we constructed 
in this work are a particular useful tool to help designing 
and testing data analysis methods for uncovering this sig- 
nal from real data. Also, we have here been able to show 
that the energy dependence of the signal offers an additional 
highly attractive way to detect specific signatures of the an- 
nihilation emission spectrum. It is quite possible that the 
different energy channels onboard the satellite FERMI that 
is currently mapping the EGB will already allow a fruitful 
application of such an approach. 
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